#include "slepbase.h"
#include "epph.h"


void nnLogisticC(slep_matrix* A, slep_vector* y, double z, OPTS* opts,
				 slep_vector* x, double *c, slep_vector* funVal){

	double sWeight[2], m1=0., m2=0., k1,k2, rsL2, lambda ; //, lambda_max;
	slep_vector *weight, *b, *x1, *Ax, *ATb, *b1;
	int i;

	assert( y->dim == A->m);

	initNormalization(A,opts);

	weight = slepAllocDenseVector(y->dim);
	b = slepAllocDenseVector(y->dim);
	x1 = slepAllocDenseVector(x->dim);
	b1 = slepAllocDenseVector(y->dim);
	ATb = slepAllocDenseVector(x->dim);
	Ax = slepAllocDenseVector(y->dim);

	if(hasGot(opts,SlepOPT_sWeight)){
		sWeight[0] = opts->sWeight[0];
		sWeight[1] = opts->sWeight[1];
		for(i=0; i < y->dim; i++)
			m1+= ((y->data.vctr[i])>0?1.:0.);
		m2 = i-m1;
		m1*=sWeight[0];
		m2*=sWeight[1];

		k1=sWeight[0]/(m1+m2);
		k2=sWeight[1]/(m1+m2);
		weight = slepAllocDenseVector(y->dim);
		m1=m2=0.0;
		for(i=0; i < weight->dim; i++){
			(weight->data.vctr[i]) = ((y->data.vctr[i])>0?k1:k2);
			m1 += ((y->data.vctr[i])>0?(weight->data.vctr[i]):0.0);
			m2 += ((y->data.vctr[i])>0?0.0:(weight->data.vctr[i]));
		}
	}else{
		weight = slepAllocDenseVector(y->dim);
		m1= 1.0 / ((double)(y->dim));
		for(i=0; i < y->dim; i++){
			m2+= ((y->data.vctr[i])>0?0.:1.);
			(weight->data.vctr[i]) = m1;
		}
		m1 = i-m2;
		m1 *= (1.0 / ((double)(y->dim)));
		m2 *= (1.0 / ((double)(y->dim)));
	}

	if( hasGot(opts,SlepOPT_rsL2)){
		rsL2 = opts->rsL2;
		assert(rsL2 > 0.);
	}else 
		rsL2 = 0.;

	//if( opts->rFlag == 0)
	//	lambda = z;
	//else{
	//	assert ( !(z<0. || z>1.) );
	//	for(i=0; i < y->dim; i++)
	//		b->data.vctr[i] = ((y->data.vctr[i])>0 ? m2 : -m1);
	//	slep_xdiagy(weight,b);

	//	NORMALIZEDMTV(A,b,ATb,opts, b1);
	//	lambda_max = slep_damax(ATb);
	//	lambda = z * lambda_max;
	//	rsL2 *= lambda_max;
	//}

	if( opts->init == 2){
		slep_zerov(x);
		(*c) = 0.0;
	}else{
		if(hasGot(opts,SlepOPT_x0))
			slep_dvcopy(opts->x0, x);
		else
			slep_zerov(x);
		if(hasGot(opts, SlepOPT_c0))
			(*c) = ((opts->c0)->data.vctr[0]);
		else
			(*c) = 0.0;
	}

	NORMALIZEDMV(A,x,Ax,opts,x1,i);

	/* The main program */ 
	{
		int bFlag=0, zf_step;
		double gamma = 1.0l/(A->m) + rsL2, lambda0;
		double cp, ccp, tp, t, alpha, sc;
		slep_vector *weighty, *Axp, *xp, *xxp;
		slep_vector *s, *aa, *As, *bb, *prob, *g, *v, *vv, *xx;
		int iterStep;
		double fun_s, gc, l_sum, r_sum, fun_x;

		
		weighty = slepCopyToNewVector(weight);
		slep_xdiagy(y,weighty);

		xp = slepCopyToNewVector(x);
		Axp = slepCopyToNewVector(Ax);
		xxp = slepAllocDenseVector(x->dim);
		slep_zerov(xxp);
		cp = (*c);
		ccp = 0;

		tp = 0.; t = 1.;
		lambda0 = 0.0;

		s = slepAllocDenseVector(x->dim);
		As = slepAllocDenseVector(Ax->dim);
		aa = slepAllocDenseVector(As->dim);
		bb = slepAllocDenseVector(As->dim);
		prob = slepAllocDenseVector(As->dim);
		g = slepAllocDenseVector(ATb->dim);
		v = slepAllocDenseVector(x->dim);
		vv = slepAllocDenseVector(x->dim);
		xx = slepAllocDenseVector(x->dim);

		for(iterStep = 0; iterStep< opts->maxIter; iterStep ++){
			alpha = (tp - 1.)/ t;
			slep_dvcopy(x,s);
			slep_daxpy(alpha,xxp,s);
			sc = (*c) + alpha * ccp;

			slep_dvcopy(Ax,As);
			slep_daxpby(-alpha,Axp,1.0+alpha, As);
			
			slep_dvcopy(As, aa);
			slep_xdiagy(y,aa);
			slep_daxpby(-sc,y,-1.0,aa);

			gc = 0.;
			fun_s = 0.;
			for(i=0; i< aa->dim; i++){
				bb->data.vctr[i] = max(aa->data.vctr[i], 0.);
				fun_s +=  weight->data.vctr[i] * ( log( 
														exp(-(bb->data.vctr[i]))+
														exp((aa->data.vctr[i])-(bb->data.vctr[i])))
													+ (bb->data.vctr[i])) ;
				prob->data.vctr[i] = 1. / (1. + exp(aa->data.vctr[i]) );
				b->data.vctr[i] = -(weighty->data.vctr[i]) *(1 - prob->data.vctr[i]);

				gc += b->data.vctr[i];
			}

			fun_s += slep_ddot(s,s) * rsL2/2.0;
			NORMALIZEDMTV(A,b,g,opts,b1);

			slep_daxpy(rsL2,s,g);

			slep_dvcopy(x,xp);
			slep_dvcopy(Ax,Axp);
			cp = (*c);

			for(;;){
				slep_dvcopy(s,v);
				slep_daxpy(-1.0/gamma, g, v);

				(*c) = sc - gc/gamma;

				//void eplb(double * x, double *root, int * steps, double * v,int n, double z, double lambda0);

				for(i=0; i< v->dim ; i++){
					vv->data.vctr[i] = (v->data.vctr[i]>=0.0? v->data.vctr[i]: 0.0);
					x->data.vctr[i] = (v->data.vctr[i]>=0.0? 1.0 : 0.0);
				}

				eplb(xx->data.vctr, &lambda, &zf_step, vv->data.vctr, vv->dim, z, lambda0);
				lambda0 = lambda;

				for(i=0; i< x->dim ; i++)
					x->data.vctr[i] = (x->data.vctr[i]>=0.0? xx->data.vctr[i] : 0.0);

				for(i=0; i < x->dim; i++)
					v->data.vctr[i] = x->data.vctr[i] - s->data.vctr[i];

				NORMALIZEDMV(A,x,Ax,opts,x1,i);

				slep_dvcopy(Ax, aa);
				slep_xdiagy(y,aa);
				slep_daxpby(-(*c),y,-1.0,aa);

				fun_x = 0.;
				for(i=0; i< aa->dim; i++){
					bb->data.vctr[i] = max(aa->data.vctr[i], 0.);
					fun_x +=  weight->data.vctr[i] * ( log( 
															exp(-(bb->data.vctr[i]))+
															exp((aa->data.vctr[i])-(bb->data.vctr[i])))
														+ (bb->data.vctr[i])) ;
				}
				fun_x += slep_ddot(x,x) * rsL2/2.0;

				r_sum = (slep_ddot(v,v) + ((*c)-sc)*((*c)-sc))/2.0;
				l_sum = fun_x - fun_s - slep_ddot(v,g) - ((*c)-sc)* gc ;

				if( r_sum <= 1e-20 ){
					bFlag=1;
					break;
				}

				if( l_sum <= r_sum * gamma)
					break;
				else
					gamma = max(2.0*gamma, l_sum/r_sum);
			}

			tp = t;
			t = (1.0 +sqrt(4.0*t*t + 1.0))/2.0;

			slep_dvcopy(x,xxp);
			slep_daxpy(-1.0, xp, xxp);
			ccp = (*c) - cp;
			funVal->data.vctr[iterStep] = fun_x ;
			if(bFlag)
				break;

			if( opts->tFlag == 0){
				if( iterStep >= 2){
					if( fabs(funVal->data.vctr[iterStep] - funVal->data.vctr[iterStep-1]) < opts->tol)
						break;
				}
			}else if( opts->tFlag == 1){
				if(iterStep>=2){
					if( fabs(funVal->data.vctr[iterStep] - funVal->data.vctr[iterStep-1]) < (opts->tol)*(funVal->data.vctr[iterStep-1]) )
						break;
				}
			}else if(opts->tFlag == 2){
				if( funVal->data.vctr[iterStep] <= opts->tol)
					break;
			}else if(opts->tFlag == 3){
				if(slep_dnrm2(xxp) <= opts->tol)
					break;
			}else if(opts->tFlag == 4){
				if( slep_dnrm2(xxp) <= (opts->tol) * max(slep_dnrm2(xp),1))
					break;
			}else if(opts->tFlag == 5){
				if( iterStep >= opts->maxIter)
					break;
			}
		}
		slepFreeVector(xx);
		slepFreeVector(vv);
		slepFreeVector(weighty);
		slepFreeVector(Axp);
		slepFreeVector(xp);
		slepFreeVector(xxp);
		slepFreeVector(s);
		slepFreeVector(aa);
		slepFreeVector(As);
		slepFreeVector(bb);
		slepFreeVector(prob);
		slepFreeVector(g);
		slepFreeVector(v);
		funVal->dim = iterStep;
	}

	slepFreeVector(weight);
	slepFreeVector(b);
	slepFreeVector(x1);
	slepFreeVector(Ax);
	slepFreeVector(ATb);
	slepFreeVector(b1);
	return;
}